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METHOD FOR SEISMIC MIGRATION USING 
EXPLICIT DEPTH EXTRAPOLATION OPERATORS 
WITH DYNAMICALLY VARIABLE OPERATOR LENGTH 

CROSS-REFERENCES TO RELATED APPLICATIONS Not Applicable 

FEDERALLY SPONSOR RESEARCH OR DEVELOPMENT Not Applicable 

SEQUENCE LISTING, TABLE, OR COMPUTER LISTING Not Applicable 

BACKGROUND OF THE INVENTION 

1 . Field of the Invention 

[0001] This invention relates generally to the field of geophysical prospecting. More 

particularly, the invention relates to the field of seismic data processing. Specifically, the 
invention is a method for seismic migration using explicit depth extrapolation operators with 
dynamically variable operator length. 

2. Description of the Related Art 

[0002] The use of three-dimensional (3D) seismic methods has resulted in increased 

drilling success in the oil and gas industry. However, 3D seismic methods are still 
computationally expensive. A crucial point in processing 3D seismic data is the migration step, 
both because of its 3D nature and the computational cost involved. The accuracy, stability, and 
efficiency of 3D migration are determined by the wavefield extrapolation technique employed. 
Thus, it is desirable that the algorithms employed in 3D migration be accurate, stable, and 
efficient themselves, especially when 3D pre-stack migration is considered. 
[0003] A number of the algorithms used in two-dimensional (2D) depth migration have 

not been successful in the 3D case. For example, implicit finite-difference extrapolation 
methods have the advantage of being unconditionally stable, but the disadvantage of bemg 
difficult to extend to the 3D case. The most common 3D implicit method is based on the 
operator splitting into alternating direction components. The splitting errors that these methods 
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exhibit in the 3D case are translated into non-circularly symmetric impulse responses, which 
become unacceptable for dips higher than 45». The splitting methods have errors that depend 
significantly on reflector dip and azimuth and thus have problems with reflector positioning 
errors Similarly, two-pass methods have problems handling lateral velocity variations. 
[0 0041 In contrast, explicit extrapolation methods that approximate the extrapolation 

operator as a finite-length spatial filter are easily extended to the 3D case. The difficulty with 
explicit extrapolation is that the stability condition is not automatically fulfilled. The stability 
condition is that no amplitude at any frequency will grow exponentially with depth. Stability 
must be guaranteed by careful design of the extrapolation operators. 

(0005] Three-dimensional seismic wavefields may be extrapolated in depth, one 

frequency at a time, by two-dimensional convolution with a circularly symmetric, frequency- and 
velocity-dependent operator. This depth extrapolation, performed for each frequency 
independently, lies at the heart of 3D finite difference depth migration. The computational 
efficiency of 3D depth migration depends directly on the efficiency of this depth extrapolation. 
For these techniques to yield reliable and interpretable results, the underlying wavefield 
extrapolation must propagate the waves through inhomogeneous media with a minimum of 
numerically induced distortion over a range of frequencies and angles of propagation. 
[00061 Holberg, 0., "Towards optimum one-way wave propagation", Geophysical 

Prospecting, Vol. 36, 1988, pp. 99-114, first disclosed explicit depth extrapolation with 
optimized operators. Holberg proposes a technique, solely for 2D depth migration, by 
generalizing conventional finite-difference expressions in the frequency-space domain. Tins 
technique yields optimized spatial symmetric convolutional operators, whose coefficients can be 
pre-computed before migration and made accessible in tables. The ratio between the temporal 
frequency and the local velocity is used to determine the correct operator at each grid point 
during the downward continuation, by fitting their spatial frequency response to the desired 
phase shift response over a range of frequencies and angles of propagation to control numencal 
distortion. Holberg's technique can be made to handle lateral velocity variations. However, the 
method only applies to 2D migration. 

[00071 Blacquiere, G., Debeye, H.W.J, Wapenaar, CPA., and Berkhout, A.J, "3-D table 

driven migration", Geophysical Prospecting, Vol. 37, 1989, pp. 925-958, extended the method 
disclosed in Holberg (1988) to 3D migration. The wavefield extrapolation is performed in the 



2 



PGS-03-03US 



space-ftequency domain as a space-dependen. spatial convolution with recursive Ktrchhoff 
extrapolation operas based on phase shift operators. The optimized operator are pre- 
ca.enla.ed and stored in a table for a range of wavenumbers. The ex.rapoia.ion is performed 
recursively in me space domain, so bom vertical and lateral velocity variations can be handled. 
The B.ac,uiere e, a., method is accurate, bu, is a full 3D method and hence computation* 
expensive. 

,0008] Hale, D., "3-D depth migration via McClellan transformaltons , Geophys.cs, Vol. 

56 No 11 (November, 1991), pp. 1778-1785, introduced a more efficient 3D scheme based on 
to McClellan transform, which gives numerically isotropic extrapolation operators. Th,s 
commonly referred to as the Hale-McClellan scheme. Given me coefficients of one-dimenstona 
frequency- and velocity-dependent filters similar to .hose used to accomplish 2D depth 
migration, McClellan transformations lead to an algorithm for 3D depm migration. Because the 
coefficients of two-dimensional depth exttapolation filters are never explicitiy computed or 
stored only .he coefficients of dre corresponding one-dimensional filters are required. Applymg 
the two-dimensional extrapolation filter increases only lureariy with me number of coefficients N 
in the corresponding one-dimensional filter, whereas tire cos. of convolution with a rwo- 
dimensional filler is generally proportional to ft. However, Hale's memod has numerical 
anisotropy. 

,00091 Two Soubaras publications: Soubaras, R, "Explicit 3-D migratron ustng 

equiripple polynomial expansion and Laplacian synthesis", 62- Ann., Interna.. M.g., Soc. Expl. 
Geophys., Expanded Abstracts, 1992, pp. 905-908, and Soubaras, R., "Explicit 3-D mixtion 
using equiripple polynomial expansion and Laplacian synthesis", GeopHysies, Vol. 61, No. 5 
(September - October, 1996), pp. 1386-1393, improved the Hale-McClellan scheme. The 
Soubaras (1992, 1996) melhod used an expansion in second-order differential operators mstead 
of the McClellan transform. Soubaras (1992, 1996) melhod also uses fine Remez algonthm ,o 
oesto .he coefficient for bo«h dre exttapolation operator and for ftre differential operators. Thrs 
memod avoids me numerical anisottopy to a large degree and is comparab.e in computational 
cost to the Hale-McClellan scheme. 

,0010] The Soubaras approach takes advantage of operator circular symmetry, as do the 

McClehan transformations, bu, avoids .he computation of a 2D filter approximating ftre cosme of 
,he wavenumbers. The approach defines a Laplacian operator, approximates dre Laplacran by 
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(he sum of two ID filters approximating the second derivatives, and approximates the exaet 
extrapolation operator by a polynomial. The seeond derivative operators and the po.ynom.al 
expansion are both calculated by the Remez exchange algorithm. 

,00111 Soflid, A., and Amtsen, B„ "Cost effective 3D one-pass depth nngrauon , 

MM Prospect Vol. 42, 1994, pp. 755-716, makes the Souharas (.992, 1996, scheme 
more cos. effective. Sollid and Amtsen (1994) used frequency adaptive optimized denvatwe 
operators. The expansion in seeond order differential operators is used, but the expans.on 
coefficients and differentia, operators are designed using a .east squares approach rather than the 
Remez algorithm. A se, of variable lengd. second order differentia, operators for each 
wavenumber has different spectra and lengths to ensure that the resulting wave extiapolahon ,s 

both accurate and efficient. n . 

,00121 Mittet, R, 2002, "Explicit 3D depth migration with a constrained operator , m 
72nd Ann. Interna.. Mtg., Soc. Expl. Geophys., Expanded Abstracts, pp. 1.48-1151, discloses a 
constrained explicit operator method which is a modification of .he full 3D scheme discussed r„ 
Blacquiere e. al. (1989), above, to make the scheme more computationally efficient. The number 
of independent operator coefficients is constiained ,o reduce tire number of computer floaong 
point operations required, thus increasing computer efficiency. The innermost coefficients m tine 
eore area of the extrapolation operator a* computed in a standard fashion. The remaimng 
outermost coefficients in tire operator, re.a.ed to very high dip and evanescent wave P ropaga.,on, 
change only as a function of radius and are constant within radial intervals. 
,00131 Thus, a need exists for an explicit depth extrapolation method for 2D and 3D 

seismic migration that is accurate, stable, and efficient. 

BRIEF SUMMARY OF THE INVENTION 

,0014, The invention is a method for seismic migration using explicit depth extrapolation 

operators with dynamically variable operator .engtit. ExpUci, depth extrapolation operators are 
constructed with variable operator lengths depending on maximum dip angle, accuracy 
condition, and wavenumber. Operator rah.es are men constructed using the explicit depth 
extrapolation operators. In a further embodiment, depth migration is performed usmg the 
explicit depth extrapolation operators from the operator rabies. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

[0015] The invention and its advantages may be more easily understood by reference to 

the following detailed description and the attached drawings, in which: 

[0016] FIG. 1 is a flowchart illustrating the processing steps of a first embodiment of the 

method of the invention for constructing explicit depth extrapolation operators with dynamically 
variable operator length; 

[0017] FIG. 2 is a flowchart illustrating the processing steps of a second embodiment of 

the method of the invention for constructing explicit depth extrapolation operators with 
dynamically variable operator length; 

[0018] FIG. 3 is a flowchart illustrating the initial processing steps for utilizing tables of 

explicit depth extrapolation operators with dynamically variable operator length, as continued in 
FIG. 4; and 

[0019] FIG. 4 is a flowchart illustrating the final processing steps for utilizing tables of 

explicit depth extrapolation operators with dynamically variable operator length, as continued 
from FIG. 3. 

[0020] While the invention will be described in connection with its preferred 

embodiments, it will be understood that the invention is not limited to these. On the contrary, the 
invention is intended to cover all alternatives, modifications, and equivalents that may be 
included within the scope of the invention, as defined by the appended claims. 

DETAILED DESCRIPTION OF THE INVENTION 

[0021] The invention is a method for seismic migration using explicit depth extrapolation 

operators with dynamically variable operator length. The method of the invention is applicable 
to the methods disclosed in the publications described above and any methods derived from or 
similar to the methods disclosed in these publications. 

[0022] Downward wavefield extrapolation transforms a seismic wavefield P(x, y, to, z) at 

lateral location x, y and depth level z to the seismic wavefield P(x, y, to, z+Az) at depth level 
z+Az by convolving P(x, y, co, z) with an extrapolation operator W(x, y, k a (x, y, z), Az). The 
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seismic wavefield P(x, y, co, z) is in the space-frequency domain, having been transformed from 
the space-time domain, if necessary, by a temporal Fourier transform. Here, x and y are the 
horizontal spatial coordinates, typically in the in-line and cross-line directions, respectively, of 
the seismic survey that collected the data. The interval Az is a spatial sampling interval or step 
length in the vertical z-coordinate direction, where depth z is measured positive downwards. A 
local wavenumber k m (x, y, z) is defined by: 



CO 



c (x, y, z) 



(1) 



where co = 2irf\s the angular frequency for frequency /and c(x, y, z) is the local propagation 
velocity of the medium at the spatial location (x, y, z). Thus, explicit depth extrapolation can be 
expressed in the space-frequency domain by the following two-dimensional spatial convolution 
along the horizontal x- and y-coordinates: 

P(x, y,co,z + Az) = W(x, y, k a ,Az)* P(x, y, co, z) 

(2) 

= \\W(x', y\ k m , Az) P(x -x',y- y', co, z) dx' dy', 

The extrapolation given in Equation (2) is applied recursively downward for all depth levels z of 
interest. One method for varying the velocity in the vertical direction is to assign different 
velocities in each depth level z. Additionally, the velocity can vary laterally in the x- and y- 
coordinate directions at each depth level z. 

[0023] The extrapolation operator W(x, y, K(x, y, z), Az) in Equation (2) is an 

approximation of the inverse spatial Fourier transform of the exact extrapolation operator D(k x , 
k y , k w , Az) for explicit depth extrapolation. The exact extrapolation operator D(k x , k y , K, Az) is 
given in the frequency-wavenumber domain by the phase shift operator: 



D(k x , k y , k„, Az) = exp(i Az J k 2 m - k] - k) ), 



(3) 
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where k x and k, are horizontal wavenumbers in the x- and y-coordinates, respectively. Strictly 
speaking the phase shift operator is only valid for homogeneous media, that is, media with no 
velocity variation. The inverse spatial Fourier transform of the phase shift operator in Equation 
(2) yields a Rayleigh operator in the space-frequency domain, which, however, must be band 
limited and truncated to be useful. As a practical alternative, the explicit extrapolation operator 
is typically approximated in the space-frequency domain. Thus, the goal of explicit depth 
extrapolation is constructing an explicit extrapolation operator W(x, y, U*. y. 4 ^ that is 
simultaneously accurate, stable, and efficient. 

[00 24] Constructing an explicit extrapolation operator W(x, y, K(x, y, z), Az) that is 

accurate means that the Fourier transform W(k x , k y , k m , Az) approximates the exact extrapolation 
operator D(k x , k y , fc. Az) given by Equation (3) closely. Thus, the difference between the 
operators D(k x , k y , U Az) and W(k x , k y , k m , Az) should be a minimum for a given frequency co 
and velocity c(x, y, z), that is, for a given k m (x, y, z). The difference between the operators D(k x , 
k k m Az) and W(k x , k y , K, Az) is measured in some appropriate norm in the frequency- 
wavenumber domain and over an appropriate range of wavenumbers, called the propagation or 
passband region of the operator. The passband region is typically defined as all wavenumbers 
less than cutoff wavenumbers k c (d x , max ) and k c (B ymaK ) given by: 

a>-sinfl„ max 6,-sinfl (4) 

KV,™)- c(w) • C (x,y,z) 

where <W and 8 ymax are the maximum dip angles selected to be accurately migrated in the x- 
coordinate and y-coordinate directions, respectively. Note that if 6 x , max = <W, then k c (d x , max ) 
k c (6 y , m J. Thus, the explicit extrapolation operator satisfies the following accuracy condition in 
the wavenumber passband region: 



1 D(k x , k y , k m , Az)-W{k x , k y , K, Az) | =mimmum, 

for 0 < k x , < k c (0 x , ma J and 0 < k y < k c (d y , max ), (5) 



.presents an appropriate norm in the frequency-wavenumber domain. The cutoff 
wavenumbers -d MW Ermine the maximum dip angle accuracy of the explicit 



where II || rep 
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extrapolation operator. 

[0025] Constructing an explicit extrapolation operator W(x, y, kjx. y, z), Az) that is stable 

means that the absolute value of the magnitude of the wavenumber response W(k x , k y> k m , Az) is 
as close as possible to unity over the wavenumber passband region and is either contained or 
suppressed for angles of propagation higher than the maximum design dip angle and in the 
evanescent region. This latter range of wavenumbers is called the stopband region of the 
operator. The stopband region is usually defined as all wavenumbers greater than or equal to the 
cutoff wavenumber and less than or equal to the Nyquist wavenumbers k x , Nyq , k y . Nyq given by: 



jc , n (6) 
A? y ' m Ay' 



where Ax and A, are spatial sampling intervals or step lengths in the horizontal x- and in- 
coordinate directions, respectively. Note that if Ax = Ay, then k, m = W Thus ' the ex P hcit 
extrapolation operator satisfies the following stability condition in the wavenumber passband 



region: 



W(k x ,k y ,k a ,Az)\ = \, 

for 0<k x ,< k c (9 x , max ) and 0 < < k c (6 y , max ), (7) 



and the following stability condition in the wavenumber stopband region: 
W(k x ,k y ,k„,Az) |<1, 



x > y ' 

for k c (6 x , max ) <k x < k x , Nyq and k c (d y , m J <k y < k yNyq (8) 

It is important that the response is stable, since the explicit extrapolation operator will be applied 
recursively. 

[0026] Constructing an explicit extrapolation operator W(x, y, k m (x, y. z), Az) that is 

efficient means that applying the operator W in a migration scheme is computationally 
inexpensive. For computational efficiency, a discrete version of Equation (2) is needed. Let 



8 



PGS-03-03US 



/, and m be integers and let Ax, Ay, and Az be step lengths in the x-, y- and z-coordinate 
directions, respectively, as defined above. Then, lateral locations can be defined discretely by x, 
= /Ax andy, =fAy. Discrete downward extrapolation transforms a discrete representation of the 
wavefield P(x b y } , co, z) at lateral location x„ yj and depth level z to the wavefield P( Xi , y h co, 
z+Az) at depth level z+Az, by convolving with a discrete extrapolation operator W(x u y } , K(x it y h 
z), Az). The discrete wavenumber k m (xu yj, *) * then given by replacing Equation (1) by the 
following discrete version: 

[0027] Analogous to the second equality in Equation (2), a discrete version of explicit 

depth extrapolation can be expressed in three dimensions by 



P(x i ,y,,^z + Az) = 

X f J W(x l ,y n ,k a> (x i ,y j ,z),te)P{x i -x l ,y j -y m ,a>,z), 

l=-L m=-M 

where L and M are called the half-lengths of the operator W in the x- and y-directions, 
respectively. The integers L and M are called operator half-lengths since the number of 
coefficients in the extrapolation operator W is (2Z+1)(2M+1). For reference, the 3D version of 
Equation (10) reduces in a straightforward manner to a 2D version given by : 



P( Xi , co, z + Az) = X *.<*i. z >' P(x ' _X " * Z) • ° 1} 



[0028] Note that the operator half-lengths, L and M, in Equations (10) and (1 1) determine 

the operator length, or size, of the explicit extrapolation operator W, and thus the efficiency of 
the extrapolation method. In principle, the operator half-lengths would be infinite for exactness. 
In practice, these operator half-lengths must be finite for computer implementation. In 
particular, these operator half-lengths should be minimized for increased computer efficiency. 
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As will be described below with reference to the present invention, the operator half-lengths can 
be made to depend optimally on the maximum wavenumber, k^fa, y jt z), at a given depth level z 
or in a sub-domain D within that depth level z. 

[0029] Holberg (1988) determines the coefficients of the transformed extrapolation 

operator W for the 2D case corresponding to Equation (1 1). However, this method is only 2D. 
Blacquiere et al. (1989) extends the method to the 3D case corresponding to Equation (10). 
However, this method is still computationally expensive. Hale (1991) introduces a more 
efficient scheme based on McClellan transforms to calculate the coefficients of the extrapolation 
operator W. However, this method has numerical anisotropy. Soubaras (1992, 1996) introduces 
an improved method using second-order differential operators instead of McClellan transforms. 
However, this method still has numerical anisotropy. Sollid, A., and Arntsen, B. (1994) use 
frequency adaptive optimized derivative operators to improve the computational efficiency of the 
Soubaras method. 

[0030] The methods for applying an extrapolation operator according to the Hale- 

McClellan approach, as discussed above for Hale (1991), Soubaras (1992, 1996), and Sollid, A., 
and Arntsen, B. (1994), are two stage algorithms. First, recursive 2D Chebyshev filters are 
applied to the wave field P(x it yj, co, z) at depth level z, resulting in auxiliary fields h(x u y Jt co, z). 
Second, the pressure at the downward extrapolated depth z+Az is expressed by: 



P0t ( .,y,,<y,z + Az) = X W l (k a (x i ,y j ,z),Az)h,(x i ,y j ,o),z) (12) 

7=0 

where W,(k m (Xi, y jt z), AzJ is an extrapolation operator now expressed in cylindrical coordinates. 
For all the above mentioned types of schemes, the extrapolation coefficients must be designed 
for a sequence of wavenumber k m values, which range from 0 to the Nyquist wavenumber, which 
is given by the maximum Nyquist wavenumber from Equation (6) or given by n, if in normalized 
coordinates. The present invention shows for these types of extrapolation schemes that the 
operator half-length L can be made to depend optimally on the maximum wavenumber k a , at a 
given depth level z or in a sub-domain D within that depth level z. This optimization yields 
increased computational efficiency. This computational efficiency is further increased with the 
constrained operator approach discussed above for Mittet, R., (2002). 
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[0031] A standard procedure is to set the operator half-length to a constant for all 

wavenumbers fe. A commonly recognized fact is that when the maximum design angle is 
increased, the operator half-length must also be increased, if the numerical accuracy is kept 
fixed. Therefore, where a maximum dip angle of 55 degrees may require an operator half-length 
of 8, a maximum dip angle of 70 may require an operator half length of 16. Increasing the 
operator half-length with a factor of 2 in this case, for a full 3D operator, increase the numerical 
work with a factor of more than 4. For example, extrapolation operator half-lengths of 3, 5, and 
12 are often associated with 30-, 50-, and 70-degree dip angle accuracies, respectively. 
[0032] There is, however, another fact that is not recognized and exploited. For a given 

maximum dip angle, keeping numerical accuracy fixed, the required operator half length varies 
with the wavenumber fe. The dependency is such that the operator half length must be increased 
with increased k*. Thus, the numerical workload of depth extrapolation can be significantly 
reduced by constructing operator tables with variable operator half lengths. 
[0033] There are several ways to obtain such operator tables. A first way to obtain the 

tables is directly in the operator optimization software. For each wavenumber k w , the strategy is 
to start with too low a value of the half-lengths L(k m ) and M(k m ), for the 3D case, or just the half- 
length L(kJ, for the 2D case. Then the value of each half-length is increased with steps of 1 
until proper convergence is achieved. This strategy is illustrated in FIG. 1, which shows a 
flowchart illustrating the processing steps of a first embodiment of the method of the invention 
for constructing tables of explicit depth extrapolation operators with dynamically variable 
operator length. 

[0034] At step 101, a maximum dip angle to be accurately migrated is selected. 

Alternatively, maximum dip angles to be accurately migrated are selected in both the x- 
coordinate and y-coordinate directions. Typically, the x-coordinate and y-coordinate directions 
are the in-line and cross-line directions, respectively, of a seismic survey used to collect seismic 
data. Alternatively, maximum dip angles in all directions are selected. 

[0035] At step 102, a type of explicit extrapolation operator for depth migration is 

selected. By way of example, but not of limitation, the method of the invention is applicable to 
the explicit extrapolation operators disclosed in the publications described above and any 
operators derived from or similar to the operators disclosed in these publications. Thus, for 
example, 3D explicit depth extrapolation can be implemented with Equation (10) and 2D explicit 



11 



PGS-03-03US 



depth extrapolation can be implemented with Equation (1 1). 

[0036] At step 103, accuracy conditions are selected for the maximum dip angles 

selected in step 101 and the type of explicit depth extrapolation operator selected in step 102. 
The accuracy conditions are selected to determine if the operator half-lengths, which determine 
the size of the explicit depth extrapolation operators, are sufficiently large. For example, the 
accuracy conditions could include the requirement that the explicit extrapolation operator 
selected satisfy the accuracy condition in the wavenumber passband region given in Equation 
(5), above. Including this accuracy condition would require calculating cutoff wavenumbers 

k Je xmax ) and MW. as § iven ^ E< * uation (4) ' t0 defme the PaSSband regi0n " CUt ° ff 
wavenumbers depend upon the maximum dip angles 6 x , max and in the x-coordinate and y- 

coordinate directions, respectively, selected in step 101. 

[0037] At step 104, a wavenumber k m is selected. The selection of the wavenumber k a is 

preferably done in a systematic manner, for computational efficiency, although systematic 
selection of the wavenumber K is not a requirement of the invention. For example, a range of 
wavenumbers k m can be taken as going from a lowest wavenumber of interest, such as zero, to a 
highest wavenumber of interest, such as a Nyquist wavenumber. Then, the selection of the 
wavenumbers k w can start with the lowest wavenumber of interest and proceed sequentially to 
the highest wavenumber of interest, or, inversely, proceed sequentially from the highest 
wavenumber of interest to the lowest wavenumber of interest. 

[0038] At step 105, an operator length is selected for the wavenumber K selected in step 

104. The operator length is preferably selected in a systematic manner, for computational 
efficiency, although systematic selection of the operator lengths is not a requirement of the 
invention. For example, the selection of the operator lengths can start with a low operator half- 
length and proceed sequentially upward by 1. 

[0039] For the 2D case, as illustrated in Equation (11), the operator length is determined 

by operator half-length L. Thus, an operator half-length L(kJ in the x-direction is preferably 
selected for the wavenumber *, selected in step 105. For the 3D case, as illustrated in Equation 
(10), the total operator length is determined by both operator half-lengths L and M. Thus an 
operator half-length L(kJ in the x-direction and an operator half-length M(kJ in the y-direction 
are preferably selected for the wavenumber K selected in step 105. 

[0040] At step 106, it is determined if the type of operator selected in step 102 with the 
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operator length selected in step 105 satisfies the accuracy conditions selected in step 103 for the 
maximum dip angles selected in step 101. For the 2D case, the operator length is preferably 
given by an operator half-length L(k m ) and for the 3D case, the operator length is preferably 
given by the pair of operator half-lengths L(kJ and M(k m ). If the answer is no, the accuracy 
conditions are not satisfied, then the process returns to step 106 to select another operator length. 
If the answer is yes, the accuracy conditions are satisfied, then the process continues to step 108. 
[0041] At step 107, it is determined if the operator length selected in step 105 is the 

smallest operator length satisfying the accuracy conditions for the wavenumber k m selected in 
step 104. For the 2D case, the operator length is preferably given by an operator half-length 
L(k m ) and for the 3D case, the operator length is preferably given by the pair of operator half- 
lengths L(k m ) and M(k a ). If the answer is no, the operator length is not the smallest satisfying the 
accuracy conditions, then the process returns to step 105 to select another operator length. If the 
answer is yes, the operator length is the smallest satisfying the accuracy conditions, then the 
process continues to step 108. 

[0042] At step 108, the operator length determined in step 105 is placed into an operator 

table. For the 2D case, the operator length is preferably given by an operator half-length L(k m ) 
and for the 3D case, the operator length is preferably given by the pair of operator half-lengths 
L(ka) and M(k m ). 

[0043] At step 109 it is determined if any wavenumbers K of interest remain to be 

selected in step 104. If the answer is yes, wavenumbers k w of interest remain, then the process 
returns to step 104 to select another wavenumber k m . If the answer is no, no wavenumbers k m of 
interest remain, then the process continues to step 1 10. 

[0044] At step 1 10, the process ends. A table of explicit depth extrapolation operators 

with dynamically variable operator length has been constructed. 

[0045] A second strategy for obtaining an operator table starts with designing a suite of 

operator tables with half-lengths ranging from minimum to maximum values of the half-lengths 
LfkJ or M(k m ), that have an error that is as small as possible for each wavenumber k w value. 
Then, for each k a value; scan the wavenumber response of this suite of operator tables with a 
fixed error criterion and accept the operator that has the shortest half-length with an acceptable 
error. In this way, a single new operator table with variable operator half-lengths L(kJ and 
M(ka>) is synthesized from a suite of fixed length operator tables. This strategy is illustrated in 
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FIG. 2, which shows a flowchart illustrating the processing steps of a second embodiment of the 
method of the invention for constructing explicit depth extrapolation operators with dynamically 
variable operator length. 

[0046] At step 201, a maximum dip angle to be accurately migrated is selected. 

Alternatively, maximum dip angles to be accurately migrated are selected in both the x- 
coordinate and y-coordinate directions. Typically, the x-coordinate and y-coordinate directions 
are the in-line and cross-line directions, respectively, of a seismic survey used to collect seismic 
data. Alternatively, maximum dip angles in all directions are selected. 

[0047] At step 202, a type of explicit extrapolation operator for depth migration is 

selected. By way of example, but not of limitation, the method of the invention is applicable to 
the explicit extrapolation operators disclosed in the publications described above and any 
operators derived from or similar to the operators disclosed in these publications. 
[0048] At step 203, accuracy conditions are selected for the maximum dip angles 

selected in step 201 and the type of explicit depth extrapolation operator selected in step 202. 
The accuracy conditions are selected to determine if the operator half-lengths, which determine 
the size of the explicit depth extrapolation operators, are sufficiently large. The accuracy 
conditions are preferably selected as described with reference to step 103 of the flowchart in 
FIG. 1. 

[0049] At step 204, preliminary operator tables are constructed for the operator type 

selected in step 202 with variable operator lengths ranging from possible minimum to maximum 
values. The preliminary operator tables are designed to have an error that is as small as possible 
for each wavenumber *. value. Thus, the accuracy conditions selected in step 203 yield a 
minimum error for each wavenumber k m . 

[0050] For the 2D case, as illustrated in Equation (1 1), the operator length is determined 

by operator half-length L. Thus, the preliminary operator tables are preferably constructed for 
operator half-lengths LM, W*J in the x-coordinate direction with smallest possible 
error for each k w . For the 3D case, as illustrated in Equation (10), the operator length is 
determined by both operator half-lengths L and M. Thus, the preliminary operator tables are 
preferably constructed for operator half-lengths L^fkJ, L max (kJ in the x-direction and 
operator half-lengths M miH (kJ, M max (kJ in the y-direction with smallest possible error for 
each k w . 
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[0051] At step 205, a wavenumber k m is selected. The selection of the wavenumber is 

preferably done in a systematic manner, for computational efficiency, although systematic 
selection of the wavenumber h, is not a requirement of the invention. For example, a range of 
wavenumbers h, can be taken as going from a lowest wavenumber of interest, such as zero, to a 
highest wavenumber of interest, such as a Nyquist wavenumber Then, the selection of the 
wavenumbers k m can start with the lowest wavenumber of interest and proceed sequentially to 
the highest wavenumber of interest, or, inversely, proceed sequentially from the highest 
wavenumber of interest to the lowest wavenumber of interest. 

[0052] At step 206, the smallest operator length satisfying the accuracy condition 

selected in step 203 for the wavenumber k m selected in step 205 is determined. The smallest 
operator length is determined by scanning the preliminary operator tables constructed in step 
204. The smallest operator length is preferably determined in a systematic manner, for 
computational efficiency, although systematic determination of the operator lengths is not a 
requirement of the invention. For example, the determination of the smallest operator length can 
start with the minimum operator length and proceed sequentially upward to the maximum 
operator length in the preliminary operator tables for the selected wavenumber k w . For the 2D 
case, the smallest operator length is preferably determined by proceeding sequentially upward 
through the operator half-lengths I m ,„^J, L^fkJ. For the 3D case, the smallest operator 
length is preferably determined by proceeding sequentially upward through the pairs of operator 
half-lengths L OT , n ffcJ, • • L max (k u ) and M mi „(kJ, . . ., M max (k w ). 

[0053] At step 207, the smallest operator length determined in step 205 is placed into an 

operator table. For the 2D case, the operator length is preferably given by the operator half- 
length L(kJ and for the 3D case, the operator length is preferably given by the pair of operator 
half-lengths L(k m ) and M(k w ). 

[0054] At step 208, it is determined if any wavenumbers ka> of interest remain to be 

selected in step 205. If the answer is yes, wavenumbers k m of interest remain, then the process 
returns to step 205 to select another wavenumber k m . If the answer is no, no wavenumbers K of 
interest remain, then the process continues to step 209. 

[0055] At step 209, the process ends. A table of explicit depth extrapolation operators 

with dynamically variable operator length has been constructed. 

[0056] In both embodiments as discussed with reference to FIGS. 1 and 2, an operator 
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table is constructing that has variable operator half-lengths L(kJ and M(kJ depending on 

wavenumber k m . For a fixed maximum propagation angle 6 max and a fixed numerical accuracy, 

the relation will be such that the operator half-length will increase with k m . 

[0057] Table 1 shows an example of how operator lengths L and M might vary for 

wavenumber K for a constrained extrapolation operator with maximum propagation angle of 65 

degrees. 





L = M 


0 < ku < 0.05 7i 


6 


0.05 n<k 0) < 0.10 7i 
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0.10 n<k 0) < 0.65 7t 
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0.65 7E<fco>< 0.75 7i 


9 


0.75 7i <^< 0.8271 


10 


0.82 n<k w < 0.88 7t 


12 


0.88 7t < *„, < 0.94 7t 


16 


0.94 71 < fco, < 71 


20 



Table 1 



[0058] The actual wavenumber intervals in Table 1 would typically vary, depending on 

such factors as operator type, maximum dip angles, and accuracy conditions. However, a 
general trend is that an operator with a short half-length can be used for wavenumber values up 
to 70 to 80 percent of the Nyquist wavenumber. Variable length operator tables require a scan of 
the velocity model before the start of wave extrapolation. This is not numerically costly and has 
to be performed only once, independent of how many frequencies there are to extrapolate. The 
simplest procedure is to determine the lowest propagation velocity at each depth level. Explicit 
depth extrapolation is performed for one co value at the time. For a given depth level z, when the 
angular frequency co and the lowest propagation velocity c mi „(z) are known, the highest 
wavenumber k w max (z) for this depth level z is given by: 
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(13) 

Thus, the required maximum extrapolation operator half-lengths L(U and M(kJ are known for 
this depth step Az. Here, the source of the gain in computer efficiency is apparent. First, for 
small and intermediate values of frequency o, high numerical accuracy can be retained with an 
extrapolation operator with short half-lengths L(kJ and M(kJ. Second, in high velocity zones, 
as, for example, in salt, the highest wavenumber KT(z) will be small due to the high value for 
the minimum propagation velocity c min (z). Again, high accuracy and high dip wave 
extrapolation can be performed with relatively short extrapolation operator half-lengths. 
[0059] A refinement of the method of the invention is to divide each depth plane z into 

several sub-domains D. Instead of determining the lowest propagation velocity c mi „(z) for the 
entire depth plane z, determine the lowest propagation velocity c min (D, z) for each sub-domain D 
within the depth plane z. Then, as in Equation (13), the highest wavenumber KT(D, z) for this 
sub-domain D within the depth level z is given for a angular frequency co by: 

Q) (14) 

This gives an additional gain in efficiency, since a relatively shorter operator is being used in 
those sub-domains that have an overall high propagation velocity. 

[0060] Now, if a node (x s , „, z) is in the domain D at the depth z, then 3D explicit depth 

extrapolation can now be implemented by using Equation (14) with Equation (10), yielding: 

P{x i ,y J ,(a,z + *z) = 

usg» * Y'\ {x „ yM y P z), Az) />(*,. -*„ yj -y m , *, Z ). 

Similarly, 2D explicit depth extrapolation can now be implemented by using Equation (14) with 
Equation (11), yielding: 
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«*— U>^» . . . x (16) 



P(x n co, z + Az) = £ Z) ' AZ) - *" Z) 



/--MC (*>.*)) 



Further, explicit depth extrapolation can now be implemented for the Hale-McClellan approach, 
as discussed above for Hale (1991), Soubaras (1992, 1996), and Sollid, A., and Arntsen, B. 
(1994), by using Equation (14) with Equation (12), yielding: 



/=0 



[0061] FIGS. 3 and 4 show two flowcharts illustrating the processing steps for seismic 

migration utilizing tables of explicit depth extrapolation operators with dynamically variable 
operator length. The operator tables are obtained by the method of the invention, preferably by 
one of the embodiments described with reference to the flowcharts in FIGS. 1 and 2. FIG. 3 
shows a flowchart illustrating the initial processing steps, to be continued in FIG. 4. 
[0062] At step 301, a seismic data set suitable for depth migration is selected. The 

seismic data set is preferably represented in the space-frequency domain, having been 
transformed from the space-time domain, if necessary, by a temporal Fourier transform. 
[0063] At step 302, a plurality of depths z are selected in the seismic data set selected in 

step 301 . Each depth z represents a depth interval in the seismic data set. The plurality of depths 
z are selected in a systematic manner, starting at the top of the seismic data set and proceeding 
sequentially downward, interval by interval. 

[0064] At step 303, a plurality of sub-domains D are selected at each depth z in the 

plurality of depths z selected in step 302 in the seismic data set selected in step 301. 
Consideration of sub-domains D at each depth z gives flexibility to the invention. Nonetheless, 
in an alternative embodiment of the method of the invention, sub-domains D at each depth z are 
not considered. For completeness and clarity, the process will be illustrated by the embodiment 
containing sub-domains, but the extension to the embodiment not containing sub-domains is 
straightforward. 

[0065] At step 304, a velocity model is determined for the seismic data set selected in 
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step 301. The velocity model determines the local velocities c(x, y, z) for the plurality of depths 
z selected in step 302 and the sub-domains D selected in step 303. The velocity model may be 
determined by any means known in the art. 

[0066] At step 305, a lowest velocity c min (D, z) is determined for each sub-domain D 

selected in step 303 at each depth z. The determination of the lowest velocity c min (D, z) in each 
sub-domain D is preferably done by scanning the velocity model determined in step 303. In an 
alternative embodiment of the method of the invention, sub-domains D at each depth z are not 
considered. A lowest velocity c min (z) is determined only at each depth z in the plurality of depths 
z selected in step 302, by scanning the velocity model determined in step 303. For completeness 
and clarity, the process will be illustrated by the embodiment containing sub-domains, but the 
extension to the embodiment not containing sub-domains is straightforward. 
[0067] At step 306, a table of explicit depth extrapolation operators with dynamically 

variable operator lengths is selected for the seismic data set selected in step 301. The table is 
preferably constructed by the method of the invention, particularly by one of the embodiments 
described with reference to the flowcharts in FIGS. 1 and 2. For the 2D case, the operator length 
is preferably given by an operator half-length L(k m ) and for the 3D case, the operator length is 
preferably given by a pair of operator half-lengths L(k m ) and M(k w ). 

[0068] The table of explicit extrapolation operators may be designed by any means 

known in the art. In particular, the method of the invention is applicable to the methods 
disclosed in the publications described above and any methods derived from or similar to the 
methods disclosed in those publications. 

[0069] At step 307, the table of explicit depth extrapolation operators with dynamically 

variable operator lengths selected in step 306 is optionally interpolated, if necessary, for the 
seismic data set selected in step 301. Interpolation allows reasonably sized operator tables to be 
constructed, stored, and reused. The interpolation can then be designed for the particular seismic 
data set velocity model used and accuracy conditions desired. 

[0070] At step 308, the process returns from step 410 in FIG. 4, if necessary, to select 

another frequency co. 

[0071] At step 309, a frequency co is selected for the seismic data set selected in step 301 . 

The selection of the frequencies co is preferably done in a systematic manner, for computational 
efficiency, although systematic selection of the frequencies co is not a requirement of the 
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invention. For example, the selection of the frequencies co can start with the lowest frequency of 
interest and proceed sequentially to the highest frequency of interest, or, inversely, proceed 
sequentially from the highest frequency of interest to the lowest frequency of interest. 
[0072] At step 310, the process returns from step 408 in FIG. 4, if necessary, to select 

another depth z. 

[0073] At step 311, a depth z is selected from the plurality of depths z selected in step 

302. The selection of the depth z is done in a systematic manner, appropriate for downward 
depth extrapolation. Thus, the depth z is selected starting at the top of the seismic data set 
selected in step 301 and proceeding sequentially downward. 

[0074] At step 312, the process returns from step 406 in FIG. 4, if necessary, to select 

another sub-domain D. 

[0075] At step 313, a sub-domain D is selected at the depth z selected in step 3 1 1 from 

the sub-domains selected in step 303. The selection of the sub-domains D is preferably done in a 
systematic manner, for computational efficiency, although systematic selection of the sub- 
domains D is not a requirement of the invention. For example, the selection of the sub-domains 
D can start at one side of the depth interval represented by the depth z and proceed sequentially 
to the other side. 

[0076] At step 314, the process continues to step 401 in FIG. 4. 

[0077] FIG. 4 shows a flowchart illustrating the final processing steps for utilizing tables 

of explicit depth extrapolation operators with dynamically variable operator length, as continued 
from FIG. 3. 

[0078] At step 40 1 , the process continues from step 3 1 4 in FIG. 3 . 

[0079] At step 402, a highest wavenumber k m max (D, z) is calculated for the sub-domain D 

selected in step 313 of FIG. 3. The highest wavenumber k m max (D, z) is calculated for the 

frequency co selected in step 309 of FIG. 3 and using the lowest velocity c min (D, z) determined in 

step 305 of FIG. 3. The calculation of the highest wavenumber k a max (D, z) is preferably by 

division of the frequency co by the lowest velocity c min (D, z), as in Equation (14). 

[0080] At step 403, a maximum operator length is selected for the sub-domain selected in 

step 313 of FIG. 3. The maximum operator length is selected on the basis of the highest 

wavenumber k u max (D, z) calculated in step 402. In the 2D case, the operator length is preferably 

selected as the operator half-length L(k a ) corresponding to the highest wavenumber k m max (D, z). 



20 



PGS-03-03US 



In the 3D case, the operator length is preferably selected as the pair of operator half-lengths 
LfkoT'fD, z)) and Mfttf^fD, z)) corresponding to the highest wavenumber k m max (D, z). 
[0081] At step 404, explicit depth extrapolation is applied in the sub-domain selected in 

step 313 of FIG. 3. Explicit depth operators are selected from the table of explicit depth 
extrapolation operators with dynamically variable operator lengths selected in step 306 of FIG. 3. 
The explicit depth extrapolation operators are selected to have no more than the maximum 
operator length selected in step 403. 

[0082] At step 405, it is determined if any sub-domains D of interest remain within the 

depth z selected in step 311 in FIG. 3. If the answer is yes, a sub-domain D of interest remains, 
then the process proceeds to step 406, to select another sub-domain D. If the answer is no, no 
sub-domain D of interest remains, then the process continues to step 407, to check the depths z. 
[0083] At step 406, the process returns to step 3 12 in FIG. 3 to select another sub-domain 

D. 

[0084] At step 407, it is determined if any depths z of interest remain within the seismic 

data set selected in step 307 in FIG. 3. If the answer is yes, a depth z of interest remains, then the 

process proceeds to step 408, to select another depth z. If the answer is no, no depth z of interest 

remains, then the process continues to step 409, to check the frequencies co. 

[0085] At step 408, the process returns to step 3 10 in FIG. 3 to select another depth z. 

[0086] At step 409, it is determined if any frequencies co of interest remain for the 

domain D selected in step 312 in FIG. 3. If the answer is yes, a frequency co of interest remains, 

then the process proceeds to step 410, to select another frequency co. If the answer is no, no 

frequency co of interest remains, then the process continues to step 41 1, to end. 

[0087] At step 410, the process returns to step 308 in FIG. 3 to select another frequency 

co. 

[0088] At step 41 1 , the process ends. 

[0089] The invention is used for 2D and 3D table driven depth extrapolation with 

optimized space-frequency domain operators. The effective number of extrapolation operator 
coefficients varies dynamically as a function of the maximum ratio of frequency to propagation- 
velocity at a depth level or, alternatively, within a sub-domain of that depth level. Use of the 
invention increases computer efficiency for explicit depth extrapolation schemes, while 
maintaining accuracy and stability. 
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[0090] It should be understood that the preceding is merely a detailed description of 

specific embodiments of this invention and that numerous changes, modifications, and 
alternatives to the disclosed embodiments can be made in accordance with the disclosure here 
without departing from the scope of the invention. The preceding description, therefore, is not 
meant to limit the scope of the invention. Rather, the scope of the invention is to be determined 
only by the appended claims and their equivalents. 
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